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We propose a new mechanism for the achievment of homochirality in life without any auto- 
catalytic production process. Our model consists of a spontaneous production together with 
a recycling cross inhibition in a closed system. It is shown that although the rate equations 
for this system predict no chiral symmetry breaking, the stochastic master equation predicts 
complete homochirality. This is because the fluctuation induced by the discreteness of pop- 
ulation numbers of participating molecules plays essential roles. This fluctuation conspires 
with the recyling cross inhibition to realize the homochirality. 
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1. Introduction and Model System 

It has long been known since the discovery of Pasteur that organic molecules in life are 
homochiral, in other words, having a completely broken chiral symmetry. 1,2 The origin of 
this homochirality remains an unsolved important puzzle. 3-6 Various mechanisms for the ger- 
mination of chirality imbalance have been proposed such as different intensities of circularly 
polarized light in a primordial era, adsorption on optically active crystals, or the parity break- 
ing in the weak interaction. 

Predicted asymmetries, however, have turned out to be very small, 5 ' 7 ' 8 and therefore 
their amplification is indispensable. Frank showed theoretically that an autocatalytic reaction 
accompanying cross inhibition can lead to the amplification of enantiomeric excess (ee) and 
to the eventual homochirality in an open system. 9 Following this work, numerous studies have 
been performed on the chiral amplification and selection in various systems. 4-6 ' 10,11 

Recently, the amplification of ee (but not homochirality) was realized in experiments 
carried out by Soai et al., 12,13 and the temporal evolution of the chemical reaction was shown 
to be explained by a second-order autocatalytic reaction. 14, 15 Stimulated by these works, we 
proposed that, in addition to the nonlinear autocatalytic reaction, a recycling process induced 
by a back reaction gives rise to the complete homochirality in a closed system. 16 Subsequently 
several theoretical works related to this mechanism have been done. 17-22 
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In these studies of chiral amplification, the autocatalytic reaction plays an essential role 
either in open systems or in closed systems. 4-6, 9-11,16-22 So far, however, any autocatalytic 
reaction has not been found in the process of polymerization, relevant for the formation 
of organic molecules in life. Granting that some pertinent autocatalytic reaction may well 
be discovered in future, it seems worthwhile to explore possibilities of chirality selection in 
non-autocatalytic way. In this paper, we demonstrate that complete chirality selection or ho- 
mochirality is possible in a closed system with spontaneous production together with recycling 
cross inhibition but without autocatalytic reaction. 

Our model consists of achiral substrate molecule A and two chiral enantiomers R and S, 
which are produced by spontaneous productions 

A^R, A^S (1) 

with the same reaction rate ko. Furthermore, R and S are assumed to react back to A as 

R + S -» 2A (2) 

with a reaction rate hq. We call this reaction a recycling cross inhibition, which looks similar 
to Frank's cross inhibition but differs in that whether R and S are recycled back in the present 
model or eliminated out of the system in the Frank's open model. 9 

In §2, the rate equation approach shows that the system has a line of fixed points and 
chirality selection is impossible. Since the fixed line is neutral in stability, the system is ex- 
pected to be susceptible to the weakest perturbations such as fluctuation. The rate equation, 
however, describes only the evolution of average quantities. To include fluctuation effect, one 
has to consider stochastic aspects of the system evolution. This feature can be taken into 
account in a stocahstic master equation approach, where the system is described by a proba- 
bility distribution function. 23-25 From stochastic analysis in §3 and §4, it will be shown that 
the fluctuation drives the system to homochirality. The effect of the fluctuation is attributed 
to the discreteness of the microscopic process, the essence of which is extracted in the system 
size expansion in §5. The result is summarized in §6. 

2. Rate equation approach 

In the rate equation approach, the reaction processes, Eq.(l) and Eq.(2), are expressed as 

^ = k a- hqts, (3) 

^ = k a - fi rs, (4) 

together with 

c = a + r + s (5) 
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where a, r, s are concentrations of species A, R, S respectively and the total concentration c is 
assumed to be constant. The conservation of total concentration expressed by Eq.(5) implies 
that R and S are recycled back to A via the cross inhibition reaction (— /xors). The trajectories 
of the evolution are easily obtained as r — s =constant, as shown by lines with arrows in Fig. 
1. The final states are obtained by solving r = s = koa — fio rs = together with a + r + s = c, 
resulting the following hyperbola 

(- + —)(- + —) = + (6) 

as is shown in Fig. 1. There, the rate for the cross inhibition n$ is chosen to be very large 
compared to that for the spontaneous production ko as cfio = 5ko, in order to draw the fixed 
line of hyperbola clearly visible away from the diagonal boundary line r + s = c. We expect, 
however, that the cross inhibition, if it exists, should be a very rare process and cfio ko- 
The conclusion of this rate equation approach is that the enantiomeric excess (ee) defined as 

* = ^ (7) 
r + s 

takes any value (— 1 < <f> < 1), depending on initial conditions, thus indicating no chirality 
selection. 




0.2 0.4 0.6 0.8 

r/c 



Fig. 1. Trajectory of the concentration r and s, prescribed by the rate equations (3) and (4), with a 
fixed line of hyperbola. The parameter is set as cfi = 5fco • 



There is, however, a subtle feature such that, along the fixed line (hyperbola), the system 
is neither stable nor unstable, namely it is neutral. Fluctuations could be decisive to destruct 
this neutrality and resolve the hyperbola into a few fixed points. In fact, a similar situation 
with a neutral fixed line appears for a closed system with a nonlinear autocatalytic process. 
We have adopted a stocahstic master equation approach for this sytem, and found chiral 
symmetry breaking in a previous paper. 25 Therefore, we adopt the same approach for the 
present system in the following sections. 
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3. Master equation approach 

For the stochastic approach, the relevant system is to be described in a microscopic way. 
Namely, the system is confined in a fixed volume V, containing species A, R, S with population 
numbers N A , Nr, Ng, respectively. Total population number of all molecules is assumed to be 
constant N = cV since the system is closed, so that we have 

N = N A + N R + Ns . (8) 

Chemical reaction is assumed to be stochastic where a microscopic state specified by the 
population number X = (N A , Nr, Ns) varies according to a certain transition probabilities 
W(X; q) of a jump q = (q A , qR, qs) to another state X' = X + q. The probability P(X, t) of 
a state X at time t then evolves according to the master equation 

dP{ Q t ,t] = E w ( x - * ^ x -«.*)-£ w ^ v) p ( x > *)■ (9) 

Q Q 

The summation by q is restricted by the conservation condition Eq.(8). The transition prob- 
abilities for the present model is explicitly expressed as 

W(N A , N R , N s ; -1, +1, 0) = k N A , 

W(N A , N R , N s ; -1, 0, +1) = k N A , 

W(N A , N R , N s ; +2, -1, -1) = fiN R N s (10) 

and those with other g's vanish. As is confirmed later, the coefficient for the microscopic cross 
inhibition fi is related to the macroscopic reaction rate /io as 

Thus the master equation in concrete form is expressed as 

dP ( NA, ^ t R,Ns,t) =k (N A + 1)[P(N A + 1,N R - 1, N S , t) + P(N A + 1, Nr, N s - 1, t)} 

+ [i(N R + l)(N s + 1)P(N A — 2, Nr + 1, N s + 1, t) 

- {2k N A + fiN R N s }P(N A , N R , N s , t). (12) 

Because of the conservation condition Eq.(8) the microscopic state of the system is spec- 
ified by two independent variables Nr and Ns in a triangular region 

0<Nr,N s ,Nr + N s <N, (13) 

as shown in Fig. 2. Transition probabilities W connect neighboring states linked along square 
edges or along a diagonal, indicated in Fig. 2. The system is equivalent to the directed random 
walk model where a random walker jumps to the right (+1 in Nr direction), to upwards (+1 
in Ns direction) and to the left-down diagonal ( — 1 in both Nr and Ns directions) in the 
Nr — Ns phase space. One then notices easily that two homochiral states (Nr, Ns) = (N, 0) 
and (0, N) are special; there is only inflow but no outflow of the walker, namely the transition 
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N R + N S = N 
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Fig. 2. Random walk in a triangular region in Nr — Ns space. A walker at a site (Nr,Ns) jumps 
along square edges or along a diagonal indicated by arrows. 



probabilities from the microscopic homochiral states (Na, Nr, Ns) = (0, N, 0) and (0,0, N) 
to other states are zero, 



W(0,N,0;q) =W(0,0, N;q) =0 



(14) 



for any jump q, so that, once the system enters into this homochiral states, it remains there. 
In this sense, the homochiral states are regarded as a sink or a kind of black hole. All other 
microscopic states are connected directly or indirectly to these two homochiral states, and 
consequently the probabilities of the other states evolve into zero, namely P(Na, Nr, Ns; t = 
oo) = for non-homo chiral states. This can be demonstrated step-by-step as follows. Because 
of Eq.(14), the master equation for the homochiral state takes the form 

dP(0,N,0,t) 



at 



--koP(l,N-l,0,t). 



(15) 



In the asymptotic limit where every temporal evolution has died out dP/dt = 0, the asymp- 
totic value for this neighboring state becomes P(l, N — 1, 0, t = oo) = 0. Furthermore, asymp- 
totic limit of the master equation for a general state X becomes 

W(X - q; q)P(X - q, oo) = £ W(X; q)P(X, oo) (16) 

Q Q 

so that if P(X , oo) = 0, then all the probabilities of states X' connected directly to the state 
X with positive W 7 s must have vanishing asymptotic values P(X', oo) = 0. Since every state 
is interconnected, this completes the proof. 

This conclusion is confirmed numerically by solving the time development of the master 
equation Eq.(12) by using the Runge-Kutta method of the fourth order. 26 Starting from a 
completely achiral initial condition with = A^ and no chiral ingredients Nr = Ns = 
0, the probability distribution, started from the initial delta-peak P(Na, Nr, Ng,t = 0) = 
Sn a ,nSn r ,oSn s ,o, remains symmetric at any time, as shown in Fig. 3. There, the redundant 
variable Na = N — Nr — Ns is suppressed, and the probability distribution is shown in 
Nr and Ns phase space. The probability contour is also shown at the basement. In this 
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(a) (b) (c) 

Fig. 3. Time evolution of the probability distribution with trajectory of the probability contour at the 
basement, obtained by numerically integrating the master equation (9) for \x = fco/10. The system 
starts from P(N, 0,0) = 1 and depicted times are at (a) kot = 1, (b)fcoi = 9 and (c)fcot = 19. The 
total population number is N = 50. 



numerical calculation, the total population number is set as N = 50 due to the limitation 
of the calculationa capacity and [i = fco/10 (or c[1q = 5ko as in Fig. 1) for the sake of 
visibility of distribution. In the early stage at a time kot = 1 in Fig. 3(a), the probability 
distribution has an approximately Gaussian shape with a central peak at the racemic fixed 
point N R = Ns = (ko/n)(y/N[i/ko + 1 — 1) 14.5 in agreement with results obtained by the 
rate eqautions. At the intermediate time k$t = 9, the probability distibution spreads along the 
hyperbola which is the fixed line found in the rate equation approach, as shown in Fig. 3(b). 
In the late stage, the probability disitribution develops sharp peaks at the two end points 
corresponding to the two homochiral states (Fig. 3(c)). 

By the use of this probability distribution, an expectation value of any function 
f (N a, N N s) at a time t is easily calculated as 

(f(N A ,N R ,N s ))t= f(N A ,N R ,Ns)P(N A ,N R ,N s ,t) (17) 

N A ,N R ,N S 

In Fig. 4, the time development of the expectation value of the population number of R 
species (N R ) t is shown. In the early stage in Fig. 4(a), the average (N R ) t increases sharply 
to the racemic fixed point value 14.5 predicted by the rate equation within a time scale of 
kot ~ 1. This development is well described by the rate equations Eq.(3) and Eq.(4) together 
with Eq.(5), as is shown by a dashed curve in Fig. 4(a). The results by the rate equations 
completely agrees with the evolution of average value (N R )t until kot ~ 1. Then, the rate 
equations predict the saturations at the racemic value, but the numerical simulation of the 
master equations indicates a slow increase. The average (N R )t approaches ultimately to the 
value N/2, corresponding to the double peak profile of the final probability distribution at 
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the two homochiral states. The slow approach is expressed well by the following form 

{N R ) t = ^-Ae- t ' T . (18) 

This exponential behavior is evident in Fig. 4(b) where the logarithmic difference ln[7V/2 — 
{Nr)A is plotted versus k$t. The fitting gives values A = 10 and l/k^r = 0.0163. 




(a) 




50 100 150 200 250 300 350 400 

k t 



(b) 



Fig. 4. Time development of the average population {Nn) t and the absolute value of the ee order 
parameter \cf>\. (a) In an early stage, k$t < 1, (Nn) t agrees with the evolution by the rate equation, 
which is shown by a dashed curve, (b) Exponential relaxation in the asymptotic region. 



The enantiomeric excess (ee) corresponds to an order parameter in a phase transition in 
standard statistical mechanics. Adopting an analogous definition of an order parameter in 
numerical simulations for magnetic phase transitions, we define the absolute value of ee as 

The time development of an ee order parameter \<ft\ is also shown in Fig. 4(a) and (b). In the very 
early stage hot < 0.2 in Fig. 4(a) \(j>\ is very large because the probability has a large amplitude 
close to the edges Nr = or Ns = 0. As time evolves, the numbers of R and S molecules 
increase and the peak position of the probability distribution leaves from both edges (Fig. 
3(a)), and thus the ee value \(f>\ drops sharply. After the peak of the probability distribution 
reaches the racemic fixed point where \(f>\ ~ 0, the probability spreads along the fixed line 
(Fig. 3(b)), and thus \(f>\ value increases steadily to a final value, unity, of homochirality. By 
plotting the logarithm of the difference, ln[l — |^(t)|], as a function of kot, as in Fig. 4(b), one 
again obtains the exponential relaxation with the same exponent l/k^T pa 0.0162. 

4. Eigenvalue Analysis 

The asymptotic relaxation of the average value and the ee order parameter turned out to 
be exponential with the same characteristic time, k§T pa 60. We consider this problem in terms 
of eigenvalues of evolution matrix of the master equation. Since master equation is a linear 
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fi/ko k,/ji 

(c) (d) 

Fig. 5. A few largest eigenvalues Aj of the time-evolution matrix plotted (a) versus total population 
A'' at a fixed cross inhibition, /u/fco = 0.1, and (c) versus ji/ko at a fixed N = 50. By plotting the 
inverse of eigenvalues in (b) and (d), one finds fco/|Aj| is linear in TV and also in ko/fi, asymptoti- 
cally. 



equation for the probability distribution, the time evolution is written by using the evolution 
matrix M 



= J2(X\M\X')P(X',t), 



dt 

where matrix elements of M are related to the transition probabilities W as 

W(X';X-X') for X'^X, 



(20) 



(21) 



(X\M\X') = < 

-E Q W(X;q) for X' = X. 
By the use of eigenfunctions and eigenvalues Aj of the matrix M as 

the time development of the probability distribution is expressed as a series 

oo 

P{X,t) = Y j a l e A > t * l . 
i=l 

Since the probability distribution satisfies the conservation as ^ x P(X, t) = 1, all the eigen- 
values must be non-positive. In the previous section, we have proven that two homochiral states 
corresponds to the final states, so that there are two degenerate zero eigenvalues Ai = A2 = 



(22) 



(23) 
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with eigenstates = Sn a ,o^n r ,nSn s ,o an d ^2 = 5n a ,o^n r ,o^n s ,n, or their linear combina- 
tions. The asymptotic temporal evolution is governed then by the third largest eigenvalue 



One can calculate eigenvalues of the matrix M numerically by using the subroutine 
"dgeev" in LAPACK. A few largest eigenvalues are shown in Fig. 5 as a function of the 
total population number ./V for a fixed cross inhibition [i = 0.1/co (Figs. 5(a) and (b)), or as 
a function of fi/ko for a fixed N = 50 (Figs. 5(c) and (d)). For [i = 0.1/co and N = 50, the 
third largest eigenvalue has a value A^/h® = —0.0163, in good agreement with the exponent 
1/koT = 0.0163 obtained in the previous section from the asymptotic final relaxation of the 
average (Nn) t and the ee values |</>(i)|. 

By plotting the inverse of eigenvalues —ko/Ai as in Figs. 5(b) and (d), one notices that 
it is linear in the total population number N and the inverse strength of the cross inhibition 
ko/fi, asymptotically. Therefore, for a very large system with a small cross inhibition, it should 
take a long time of the order N/ [i before the homochirality becomes observable. 

5. System size expansion 

For our model with spontaneous production of chiral species with recycling cross inhibition, 
the rate equation tells us no chirality selection whereas the stochastic master equation insists 
the final configuration be homochiral. The totally different conclusions are ascribed to the 
fluctuation due to the discreteness of microscopic processes. The fluctuation effect associated 
to the system size is qualitatively analyzed by the system size (precisely said, the inverse 
system size) expansion analysis of the master equation, developed by R. Kubo et al. 27 ' 28 

In the master Eq.(9), the probability density P(X, t) of a microscopic state X is connected 
to another state X + q which differs with a jump q of order unity by a transition probability 
W(X; q). The rate W is of macroscopic order of the system size N or the volume V as 



Then, the probability is assumed to take the form P(X,t) = exp[Vx(x, t)], and time evo- 
lutions of the leading order contributions of the average density (x) t and the correlation 
functions 



As- 



W(X;q) = Vw(x;q) 



(24) 



where x is the density variable x = X /V of order unity: 




(25) 



<Tij(t) = V (( x i ~ ( x i)t){xj - (Xj) t )) 



(26) 



are shown to be determined by moments 




(27) 



<3 
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of the transition probability w as 

ft( x i)t = c i(( x )t) 

= E (a&^W + °*®$tu) + (28) 

The evolution equations for the higher order corrections can also be derived. 27 If the system 
is normal, the fluctuation correlations o; L j remains of order unity. On the other hand, if the 
system is unstable, as in the case of phase transitions with critical behaviors, at least one of 
the fluctuations er^- is enhanced to the order of the system size. 
In the present model, first order moments are 

c r = c s = k$a — nors = ko(c — r — s) — [i$rs (29) 

and the second order moments are 

c rr = c ss = k a + fi rs = k (c - r - s) + fi rs, c rs = [i rs, (30) 

where (1q = V [i as defined previously in Eq.(ll). Thus, the lowest order of the average concen- 
trations (r)t and (s)t satisfy the rate equations Eq.(3) and Eq.(4). For simplicity, we describe 
these lowest order average values as r and s, hereafter. The correlation functions follow the 
evolution 

^cr rr (t) = -2(k + iiQs)a rr - 2(k + / u r)cj rs + k (c - r - s) + ^rs 
-^ a rs(t) = ~(ko + W)s)oYr - (2k + W + fM)s)a rs ~ {h + Hor)a ss + fi rs 

^. a ss(t) = -2(fc + Vor)a ss - 2(k + fios)a rs + k (c - r - s) + /J, rs (31) 

To detect the chiral symmetry breaking, it is more convenient to use the following symmetric 
and asymmetric variables 

x + = r + s, X- = r — s. (32) 
Their averages and correlation functions evolve as 

-^x+ = 2k (c - x+) - ^M x + ~ x -), = 0, 



= -2(2fc + fi x+)v++ + 2k (c - x + ) + w{x\ - x 2 _), 



d_ 
~dt 

^ a +- = - (2fco + HoX + )a + - + fx x-a — , — = 2A; (c - (33) 

By starting from the completely achiral initial condition, the system remains racemic 
with average value x_ = 0, and x + approaches to the racemic fixed point value 



(2/co//io)(\/ c / u o/^o + 1 — 1)- The time the system takes to reach the racemic fixed point is 
of order l/ko- 
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The fluctuation of the asymmetry variable a , however, increases indefinitely in this 

approximation. In fact, when the average value of x+ takes the racemic point value, the 
fluctuation increases linearly in time with a positive velocity; 

This increase of the fluctuation indicates that the discreteness in the microscopic process of 
chemical reaction evokes an intrinsic instability of the racemic fixed point (and in general, 
every points on the fixed line) to macroscopic level. 





(a) 



(b) 



Fig. 6. Fluctuation correlation functions obtained by numerically time integrating the master equa- 
tion (a) at an early stage, and (b) a whole time range. Parameter is set as [i — fc /10, or c^ = 5/co- 
A dotted line in (a) represents the result of system size expansion. 



In the numerical simulation the fluctuation of asymmetric variable a (t)/c actually in- 
creases in time, as shown in Fig. 6. The initial increase shown by a continuous curve in Fig. 
6(a) is linear in time hot, as expected from the system size expansion analysis, shown by a 
dashed line. As time passes, the system size analysis fails, since the average value x+ devi- 
ates from the racemic value in the actual simulation. One then has to consider higher order 
contributions to the average of x+. Departure from the racemic state becomes visible in a 

macrosopic sense when the fluctuation expressed by a /c reaches the order N. The time for 

this to appear is of order as is obtained from Eq.(34) and by assuming x + ~ c. 

Asymptotically, the double peak structure developes in the probability distribution. The 
time scale for this to fully develop is governed by the largest nonzero eigenvalue A3 so that the 
homochirality is realized in the time of order of 1 /| A3 1 ss N/fj,, as is described in the previous 
section. 

6. Conclusion and discussion 

It seems to be a general consensus that an autocatalytic production process, either linear 
or nonlinear, is indispensable for the realization of homochirality. In the present work, we 
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proposed a new mechanism by using a simple model to demonstrated that the homochirality 
can be realized without any autocatalytic production process. Our model consists of a spon- 
taeous production together with a recycling cross inhibition in a closed system. As was shown, 
the rate equations for this system predict no chiral symmetry breaking, but the stochastic 
master equation predicts complete homochirality. This is because the fluctuation induced by 
the discreteness of population numbers of participating molecules plays essential roles. This 
fluctuation conspires with the recyling cross inhibition to realize the homochirality. 

If this fluctuation mechanism could explain the homochirality in life, then this is what 
Pearson suggested long time ago. 29 However, the necessary time for the homochirality to 
set in due to the fluctuation is very large , as it is proportional to the total number of the 
relevant molecules (~ N/fi), which is of macroscopic size. Taking this feature into account, 
we can conceive a new senario for the homochirality in macroscopic scale as follows. At first, 
in some small closed corner, such as in a region enclosed by a vesicle, the fluctuation induced 
homochirality is realized in a very long time with respect to the time scale of laboratory 
experiments, but not so long with respect to geological time scale. It may be conceivable 
that in this period some unknown autocatalytic reaction, the effect of which is too small to 
be detected in the laboratory experiment so far, begins to operate and generates the large 
scale realization of the homochiralty from the homochiral seeds produced by the fluctuation 
mechanism proposed here. 

Even though reaction systems with a recycling cross inhibition are not yet found, we hope 
that a simple system with only a recycling cross inhibition might be found in near future, and 
the establishment of homochirality be checked. 
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